arXiv: 1509.07083v2 [quant-ph] 24 Sep 2015 


Time-dependent Hamiltonian estimation for Doppler velocimetry of trapped ions 


L. E. de Clercq, R. Oswald, C. Fliihmann, B. Keitch, D. Kienzler, 

H.-Y. Lo, M. Marinelli, D. Nadlinger, V. Negnevitsky, J. P. Home 
Institute for Quantum Electronics, ETH Zurich, Otto-Stern-Weg 1, 8093 Zurich, Switzerland 


The time evolution of a closed quantum sys¬ 
tem is connected to its Hamiltonian through 
Schrodinger’s equation. The ability to estimate 
the Hamiltonian is critical to our understand¬ 
ing of quantum systems, and allows optimization 
of control. Though spectroscopic methods al¬ 
low time-independent Hamiltonians to be recov¬ 
ered, for time-dependent Hamiltonians this task 
is more challenging mm- Here, using a single 
trapped ion, we experimentally demonstrate a 
method for estimating a time-dependent Hamil¬ 
tonian of a single qubit. The method involves 
measuring the time evolution of the qubit in a 
fixed basis as a function of a time-independent 
offset term added to the Hamiltonian. In our 
system the initially unknown Hamiltonian arises 
from transporting an ion through a static, near¬ 
resonant laser beam (3* Hamiltonian estimation 
allows us to estimate the spatial dependence of 
the laser beam intensity and the ion’s velocity as 
a function of time. This work is of direct value 
in optimizing transport operations and transport- 
based gates in scalable trapped ion quantum in¬ 
formation processing [S HlOj . while the estimation 
technique is general enough that it can be ap¬ 
plied to other quantum systems, aiding the pur¬ 
suit of high operational fidelities in quantum con¬ 
trol [mug. 

Estimation of the underlying dynamics which drive the 
evolution of systems is a key problem in many areas of 
physics and engineering. This knowledge allows control 
inputs to be designed which account for imperfections 
in the physical implementation. Eor closed quantum 
systems, the time dependence of a system is driven by 
the Hamiltonian through Schrodinger’s equation. If the 
Hamiltonian is static in time, a wide range of techniques 
have been proposed for recovering the Hamiltonian 
ElIIH], which have been applied to a variety of systems 
including chemical processes [4] and quantum dots [SlISl. 
These methods often involve estimation of the eigenvec¬ 
tors and eigenvalues of the Hamiltonian via spectroscopy, 
or through pulse-probe techniques for which a Fourier 
transform of the time-evolution gives information about 
the spectrum. However these methods are not directly 
applicable to time-dependent Hamiltonians. Such Hamil¬ 
tonians are becoming of increasingly important as quan¬ 
tum engineering pursues a combination of high opera¬ 
tional fidelities and speed, often involving fast variation 
of control fields which are particularly susceptible to dis¬ 
tortion before reaching the quantum device [siiniiiiHii]. 

In this Letter, we propose and demonstrate a method 


for reconstructing a general time-dependent Hamiltonian 
with two non-commuting terms which drives the evolu¬ 
tion of a single qubifi The method works with any single 
qubit Hamiltonian H = where the fi{t) are 

arbitrary time-dependent functions and are the Pauli 
operators. In our experiments, a Hamiltonian with two 
non-commuting time-dependent terms arises when we try 
to perform quantum logic gates by transporting an ion 
through a static laser beam m ng. In this case, the 
Hamiltonian describing the interaction between the ion 
and the laser can be written in an appropriate rotating 
frame as 

Hl{t) = ^ (-n(i)(Ta; + S{t)az ) (1) 

which includes a time-varying Rabi frequency and 

an effective detuning S{t) which is related to the first- 
order Doppler shift of the laser in the rest frame of the 
moving ion (see [19] for details). For a Hamiltonian of 
this type with unspecified time-dependent coefficients, no 
analytical solution to Schrodinger’s equation exists m 
\2T\ . In order to reconstruct the Hamiltonian we make use 
of two additional features of our experiment. The first is 
that we can switch off the Hamiltonian at time toff on a 
timescale which is fast compared to the evolution of the 
qubit. Secondly we are able to offset one of the terms in 
the Hamiltonian, in our case by adding a static detuning 
term Hg = HSl^z I such that the total Hamiltonian is 
Hi{t) -|- Hg. We then measure the expectation value of 
the qubit in the az basis as a function of 6 l and toff. 
Repeating the experiment with identical settings many 
times, we obtain an estimate of the expectation value 
which we denote as (^off:^L))- 

Hamiltonian extraction involves theoretically generat¬ 
ing the qubit populations (d™ {toS^^L))^ and attempt¬ 
ing to find the Hamiltonian for which this most closely 
matches the data. In order to provide a simple parame¬ 
terization, we represent 6{t) and fl{t) as a linear weighted 
combination of basis splines [22][23|. (d|^™ (toff, ^l)) is 
compared to the measured data using a weighted least- 
squares cost function, which we optimize with respect 
to the weights of the basis-splines used to parameterize 
S{t) and 0(t). Solving this optimization problem in gen¬ 
eral is hard because the cost function is subject to strong 
constraints imposed by quantum mechanics, producing 
a non-trivial relation between the weights and the spin 
populations m- We overcome this problem by making 
use of the inherent causality of the quantum-mechanical 
evolution, and by assuming that the parameters of the 
Hamiltonian vary smoothly. We call our technique “Ex¬ 
tending the Horizon Estimation”, in analogy to estab- 
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lished methods in engineering [24] (a detailed description 
of our method can be found in USD- Rather than opti¬ 
mizing over the whole data set at once, we build up the 
solution by initially fitting the data over a limited region 
of time 0 < toff < Tq. The solution obtained over this 
first region can be extrapolated over a larger time span 
0 < t < Ti where Ti = Tq +r, which we use as a starting 
point to find an optimal solution for this extended region. 
This procedure is iterated until = max(toff). The 

method allows us to choose a reduced number of basis 
spline functions to represent S{t) and r^(t), and also re¬ 
duces the amount of data considered in the early stages 
of the fit, when the least is known about the parame¬ 
ters. This facilitates the use of non-linear minimization 
routines, which are based on local linearization of the 
problem and converge faster near the optimum. More 
details regarding the optimization routine can be found 
in HD. 

In the experimental work, we demonstrate reconstruc¬ 
tion of the spin Hamiltonian for an ion transported 
through a near-resonant laser beam. Our qubit is en¬ 
coded in the electronic states of a trapped calcium 
ion, which is defined by |0) = \‘^Si/ 2 ,Mj = 1/2) and 
|1) = |^T) 5 / 2 , Mj = 3/2). This transition is well resolved 
from all other transitions, and has an optical frequency 
uJo/{27r) 411.0420 THz. The laser beam points at 

45 degrees to the transport axis, and has an approxi¬ 
mately Gaussian spatial intensity distribution. The time- 
dependent velocity z{t) of the ion is controlled by adia¬ 
batic translation of the potential well in which the ion is 
trapped. This is implemented by applying time-varying 
potentials to multiple electrodes of a segmented ion 
trap, which are generated using a multi-channel arbitrary 
waveform generator, each output of which is connected 
to a pair of electrodes via a passive third order low-pass 
Butterworth filter. The result is that the ion experiences 
a time-varying Rabi frequency Q{t) and a laser phase 
which varies with time as 4>(t) = (j){z{t)) — where 
(j){z{t)) = kz{z{t))z{t) with kz{z{t)) the laser wavevector 
projected onto the transport axis at position z{t) and ujl 
the laser frequency. The spatial variation of kz{z{t)) ac¬ 
counts for the curvature of the wavefronts of the Gaussian 
laser beam. In order to create a Hamiltonian of the form 
of equation we work with the differential of the phase, 
which gives a detuning S{t) = Sl — 4> = {k'^{z)z -h kz{z)) z 
with 5l = ujl—^o the laser detuning from resonance. For 
planar wavefronts k'^{z) = 0, and 6{t) corresponds to the 
familiar expression for the first-order Doppler shift (see 
m for details). 

The experimental sequence is depicted in figure We 
start by cooling all motional modes of the ion to n < 3 
using a combination of Doppler and electromagnetically- 
induced-transparency cooling [25| , and then initialize the 
internal state by optical pumping into |0). The ion is 
then transported to zone A, and the laser beam used to 
implement the Hamiltonian is turned on in zone B. The 
ion is then transported through this laser beam to zone C. 
During the passage through the laser beam, we rapidly 




FIG. 1: Experimental sequence and timing: a) The ex¬ 
periment is carried out in three zones of the trap indicated 
by A, B and C. b) The experimental sequence involves steps 
(i) through (v). Preparation and readout are carried out on 
the static ion in zone B. The qubit evolves while the ion is 
transported from zone A to zone C, via the laser beam in 
zone B. c) Experimental sequence showing the timing of ap¬ 
plied laser beams and ion transport, including shutting off the 
laser beam during transport. 


turn the beam off at time toff and thus stop the qubit 
dynamics. The ion is then returned to the central zone 
B in order to perform state readout, which measures the 
qubit in the computational basis (for more details see 
USD- The additional Hamiltonian is implemented 
by offsetting the laser frequency used in the experiment 
by a detuning Sl- For each setting of toff and Sl the 
experiment is repeated 100 times, allowing us to obtain 
an estimate for the qubit populations (toff,<^L))- 
Experimental data is shown in figurefor two different 
beam positions, alongside the results of fitting performed 
using our iterative method. The beam positions used for 
each data set differ by around 64 pm, but the transport 
waveform used was identical. The reconstructed veloci¬ 
ties should therefore agree in the region where the data 
overlap. It can be seen from the residuals that the es¬ 
timation is able to find a Hamiltonian which results in 
a close match to the data. In order to get an estimate 
of the relevant error bars for our reconstruction, we have 
performed non-parametric resampling with replacement, 
optimizing for the solution using the same set of B-spline 
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FIG. 2: Measured data, estimation and residuals: Spin population as a function of detuning and switch-off time of the 
laser beam, a) is for a laser beam centered in zone B, while for b) the beam was displaced towards zone C by 64 /xm. From 
left to right are plotted the experimental data, the populations generated from the best fit Hamiltonian, and the residuals. 
Each data point results from 100 repetitions of the experimental sequence. The data in a) consist of an array of 100 x 101 
experimental settings, while that shown in b) consists of an array of 201 x 201 settings. This leads to smaller error bars in 
the reconstructed Hamiltonian for the latter. For the Hamiltonian estimation the data was weighted according to quantum 
projection noise. 


functions as was used for the experimental data to pro¬ 
vide a new estimate for the Hamiltonian. This is repeated 
for a large number of samples, resulting in a distribution 
for the estimated values of 6{t) and from which we 
extract statistical properties such as the standard error. 
The error bounds shown in figures correspond to the 
standard error on the mean obtained from these distri¬ 
butions (see m for further details). 

The estimated coefficients of the Hamiltonian ex¬ 
tracted from the two data sets are shown in figure]^). 
It can be seen that the values of 6{t) for the two different 
beam positions differ for the region where the reconstruc¬ 
tions overlap. We think that this effect arises from the 
non-planar wavefronts of the laser beam. Inverting the 
expression for 6{t) to obtain the velocity of the ion, we 
find z{t) = 5(t)/(k'^{z)z + kz{z)). Using this correction, 
we find that the two velocity profiles agree if we assume 
that the ion passes through the center of the beam at a 
distance of 2.27 mm before the minimum beam waist, a 
value which is consistent with experimental uncertainties 
due to beam propagation and possible mis-positioning of 


the ion trap with respect to the fixed final focusing lens. 
The velocity estimates taking account of this effect are 
shown in figure]^). 

Figure shows the results of a reconstruction for a sec¬ 
ond pair of data sets taken using two different velocity 
profiles but with a common beam position. The resolu¬ 
tion in both time and detuning were lower in this case 
than for the data shown in figure]^ (see m for the data). 
We observe that the estimated Rabi frequency profiles 
agree to within the error bars of the reconstruction. One 
interesting feature of this plot is that the error bars pro¬ 
duced from the resampled data sets are notably higher at 
the peak than on the sides of the beam. We think that 
this happens because the sampling time of the data is 
0.5 /xs, which is not high enough to accurately resolve the 
fast population dynamics resulting from the high Rabi 
frequency (the Nyquist frequency is 1 MHz). In order to 
optimize the efficiency of our method, it would be advan¬ 
tageous to run the reconstruction method in parallel with 
data taking, thus allowing updating of the sampling time 
and frequency resolution of points based on the current 
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FIG. 3: Estimates of time-dependent co-efRcients: a) 

The effective detuning S{t) and c) Rabi frequency Q(t) ob¬ 
tained from the two data sets, along with dashed lines in¬ 
dicating the standard error on the mean of these estimates 
obtained from resampling. For a), the inset shows a close 
up of the estimated S{t) in the regions where the estimates 
overlap, showing that these do not give the same value, b) 
The estimated velocity z{t) of the ion obtained after applying 
wavefront correction. The inset shows that this can produce 
consistent results. 




position (/im) 

FIG. 4: Spatial Rabi frequency: a) The estimated 6{t) 
obtained from the second pair of data sets (Figure]^ in US]), 
b) The estimated Rabi frequency Q(t) for the same two data 
sets. 


estimates of parameter values. 


Our method for directly obtaining a non-commuting 
time-dependent Hamiltonian uses straightforward mea¬ 
surements of the qubit state in a fixed basis as a func¬ 
tion of time and a controlled offset to the Hamiltonian. 
This simplicity means that the method should be appli¬ 
cable in a wide range of physical systems where such con¬ 
trol is available, including many technologies considered 
for quantum computation m 0 El EH HD. A process- 
tomography based approach would require that for ev¬ 
ery time step multiple input states be introduced, and a 
measurement made in multiple bases [28U3Q] . An effec¬ 
tive modulation of the measurement basis arises in our 
approach due to the additional detuning Sl- It is worth 
noting that tomography provides more information than 
our method: it makes no assumptions about the dynam¬ 
ics aside from that of a completely positive map while 
we require coherent dynamics. Extensions to our work 
are required in order to provide a rigorous estimation of 
the efficiency of the method in terms of the precision ob¬ 
tained for a given number of measurements, and to see 
whether a similar approach could be taken to non-unitary 
dynamics. Using this method on considerably lower res¬ 
olution data sets, we have recently been able to improve 
the control over the velocity, which will be necessary in 
order to realize multi-qubit transport gates in our current 
setup [7]- 
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FIG. 5: Beam and ion transport: The beam propagation 
direction lies along the ^-axis and the ion is transported along 
the z-axis lying the /^^-plane as indicated. Normalized vectors 
representing ei{K, lying perpendicular to the wavefronts are 
indicated by the blue arrows. 


I. SUPPLEMENTARY MATERIAL 
A. Derivation of Hamiltonian 

The interaction of a laser beam with frequency cjl 
and wave vector with a two-level atom with res¬ 

onant frequency ujq and time-dependent position of the 
ion z{t) = (0, 0, z{t)) can be described in the Schrodinger 
picture by the Hamiltonian 



where the Rabi frequency ft{z{t)) gives the interaction 
strength between the laser and the two atomic levels. We 
can define the laser phase at the position of the ion as 
^(t) = 0(t) —uJLt with 0(t) = k{z{t)) • z{t) = kz{z{t))z{t) 
and kz{z{t)) = \k\cos{0{t)) being the projection of the 
laser beam onto the ^-axis along which the ion is trans¬ 
ported. Here 0{t) is the angle between the wave-vector 
k{z{t)) and the transport axis evaluated at position z{t). 
Moving to a rotating frame using the unitary transforma- 

_ ■ <^(t) 

tion U = e * 2 and applying the rotating wave approx¬ 
imation with respect to optical frequencies, we obtain 

Hi = ^ j az ) . (3) 

Defining a static detuning Sl = ojl — obtain 

Hi = ^ + (V - Oz ) . (4) 

with 

S{t) =Sl - (pit) (5) 

which is the expression used in the main text. 


B. Wavefront correction 

For plane waves we find that (p{t) = k • v{t) which 
is the well-known expression for the first-order Doppler 


shift. For transport through a real Gaussian beam, the 
wave-vector direction changes with position. Taking this 
into account, the derivative of 0(t) becomes 


Ht) = + kziz{t))] z{t) (6) 


where k'^ = dkzjdz and z{t) is the component of the ion’s 
velocity which lies along the 2 ;-axis. We extract 6{t) using 
our Hamiltonian estimation procedure, thus to obtain the 
velocity of the ion we use 


z{t) 


-S{t) + Sl 

K(z(t))z{t) + kz{z{t)) 


( 7 ) 


As the ion moves through the beam it experiences the 
same magnitude of the wave vector |^| = 27r/A, but the 
angle 0 between the ion’s direction and the wave vector 
changes. Written as a function of this angle, the velocity 
becomes 


z{t) 


_ -S{t) + Sl _ 

—1^1 sin {0{z{t))) 6'{z{t))z{t) + |^| cos(6>(^(t))) 


( 8 ) 


where 0'{z{t)) = d0{z{t))/dz{t). We parameterize our 
Gaussian beam according to figure The phase is given 
as a function of both the position along the beam axis ^ 
and the perpendicular distance from this axis hz by m 




\k\K‘^ 

m^) ■ 


(9) 


where the Gaussian beam parameters include the beam 
waist W{^), the radius of curvature R{^) the Rayleigh 
range and the Guoy phase shift C(^). These are given 
by the expressions 


W{0 = Wo^l+(^^^ 

77(0 = ^ (^1 + (y) j 

ao = ta„-‘(i-) 

ttWS 
^H = — 


( 10 ) 


where Wo is the minimum beam waist and A the laser 
wavelength. The ion moves along the z-axis shown in 
figure 1^ In the ^c^-plane a unit vector perpendic¬ 

ular to the wavefronts is given by 


ei{K,0 


iiv^(«,eii 


(11) 


and the unit vector pointing along the direction of 
transport is given by 


cos{a) 

sin(a) 


( 12 ) 
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The angle 0{^) between the wave and position vector 
is then given by the dot product 


which can be written in terms of the full set of parameters 
above as 


0{k) = cos ^ (e^ • ey). (13) 


0{k,) = cos ^ (71 + 72) 

cos(a) (-2^r + kK^ (Cl - C^) + 2k (C^ + C|)^) 


7i = 


ri{K) 


72 = 


sm{a)2kK^ (C^ + C|) 


7?(k) 


= if+^ r) 

K{t) = z{t) sin((a) 




f + en V if + fR) . 


(14) 




time 


FIG. 6: Extending the Horizon Estimation: The steps 
performed when extending the time horizon from Tn to Tn+i 
are illustrated. We first predict in the old basis, then move 
to the new basis, and finally optimize again. The figure also 
shows the basis splines 


and fl{t) with basis spline curves. Basis spline curves 
allow the construction of smooth functions using only a 
few parameters. This is achieved by introducing a set of 
polynomial Basis (B)-spline functions of order k 

[23] , A smooth curve S{t) can then be represented as a 
linear combination of these B-spline basis functions [22] 


= ( 15 ) 

i=0 


The B-splines Bi^k{t) of order k are recursively defined 
over the index i over a set of points K = {to, ti,..., tn-\-k} 
which is referred to as the knot vector [23] . 


B, 


U <t < ti+i 
otherwise 


ift) = {J 

= uJi,k{t)Bi^k-i{t) + (1 


_ 

ti+k 

0 


^ ^i+k — l 

otherwise 

(16) 


where in our experiments a = 37r/4. 

Using Eq. and we examined the value of ^ci re¬ 
quired for the velocity to match for our two beam posi¬ 
tions. We find that they agree for ^ci = —2.27 mm, which 
is within the experimental uncertainties for our setup. 


C. Basis spline curves 

One challenge in obtaining an estimate for the Hamil¬ 
tonian is that we must optimize over continuous func¬ 
tions S{t) and Q{t). To address this, we represent S{t) 


Figure]^ gives a visualization of the B-splines Bi^kif) and 
a basis spline curve. The B-spline construction ensures 
that any linear combination of the B-splines is continuous 
and has {k—2) continuous derivatives. The knot vector K 
determines how the basis functions are positioned within 
the interval We notice that for our Hamilto¬ 

nian the spacing of the B-splines is not critical, which 
we think is due to the smoothness of the variations in 
our Hamiltonian parameters S{t) and Q{t). We therefore 
used the Matlab function spap2 to automatically choose 
a suitable knot vector and restricted ourselves to opti¬ 
mizing the coefficients ai. We collect all coefficients ai 
for S{t) and Q{t) and store them in a single vector a. 























D. Extending the Horizon Estimation 


The task of inferring the time-dependent Hamiltonian 
of the form 0 from the measured data can be cast into 
an optimization problem for which we use a reduced chi- 
squared cost function 




toff 


{toS,SL)) - itoS,SL)) 






(17) 


where iy = N — n — lis the degrees of freedom with N 
the number of data points and n the number of fitting 
parameters, and is the standard error on 

the estimated (toff, Sl)) which we obtain assuming 

quantum projection noise. In our case the Hamiltonian 
Hi (Eq. Ih with offset Hg = HSlO'z I is parametrized by 
fl{t) and 6{t). We can thus write the problem as 

^min (18) 

subject to 

|«'(i = 0,(5L)) = |0), 

(<7®™ (i,(5L)) = \'f>{t,6L)) (19) 

for all 6l- 

This optimization problem is hard to efficiently solve 
in general, because it is nonlinear and non-convex due 
to the nature of Schrodinger’s equation and the use of 
projective measurements. In order to overcome this chal¬ 
lenge, we have implemented a method which we call “Ex¬ 
tending the Horizon Estimation” (EHE) in analogy to a 
well-established technique called “Moving Horizon Esti¬ 
mation” (MHE) [24] . 

The key idea is that because our measurement data 
arises from a causal evolution, we can also estimate the 
Hamiltonian in a causal way. We define a time span 
ranging from the initial time to some later time which 
we call the time horizon. Instead of optimizing J over 
the complete time span at once, we first restrict ourselves 
to a small, initial time horizon reaching only up to the 
start of the qubit dynamics. Optimizing J over this short 
time horizon requires fewer optimization parameters and 
is simpler than attempting to optimize over the full data 
set. Once we have solved this small sub-problem, we 
extend the time horizon and re-run the optimization, ex¬ 
trapolating the results of the initial time window into the 
extended window in order to provide good starting con¬ 
ditions for the subsequent optimization. This is greatly 
advantageous for the use of non-linear least squares opti¬ 
mization, which typically works by linearizing the prob¬ 
lem and converges much faster near the optimum. The 
extension of the horizon is used repeatedly until the time 
window covers the full data set. 

Conceptually EHE is very similar to MHE. The main 
difference is that in MHE the time span has a fixed length 


and thus its origin gets shifted forwards in time along 
with the horizon. In EHE the origin stays fixed at the 
expense of having to increase the time span under con¬ 
sideration. MHE avoids this by introducing a so-called 
arrival cost to approximate the previous costs incurred 
before the start of the time span. This keeps the compu¬ 
tational burden fixed over time, which is very important 
as MHE is usually used to estimate the state of a sys¬ 
tem in real-time, often on severely constrained embed¬ 
ded platforms. Since neither constraint applies to our 
problem, we decided to extend the horizon rather than 
finding an approximate arrival cost. This is advantageous 
since finding the arrival cost in the general case is still 
an open problem. Due to the similarity between MHE 
and EHE, we anticipate future improvements by adapt¬ 
ing techniques used in MHE to EHE. 

Next, we present a more detailed algorithmic summary 
of our implementation of the method outlined above. 

1. Searching for a starting point. Here we re¬ 
construct the Hamiltonian for a first, minimal time 
horizon such that we can then use this as a starting 
point to iteratively extend the horizon as described 
in step 2. 

(a) Choose an initial time horizon such that it 
contains the region where the first discernible 
qubit dynamics occur. 

(b) Cut down the number of fitting parameters 
as much as possible, e.g. by using few Basis 
splines of low order. This amounts to choosing 
empirically a low number of basis splines (and 
thus the length of ^o) which might represent 
5{t) and 0(t) over the given region. 

(c) Use a nonlinear least squares fitting routine 
to minimize J by varying the parameters ao. 
In the case that the initial fit is not good or 
no minimum is found, try new initial condi¬ 
tions, change the number of B-spline func¬ 
tions, or manually adjust the function using 
prior knowledge of the physical system under 
consideration. 

This procedure is used to provide a starting point 
for the optimization over the initially chosen win¬ 
dow, which is typically performed with a higher 
order set of B-splines. Erom this starting point, 
we iteratively extend the fitting method to the full 
data set as follows. 

2. Extend the horizon This step is repeated until 
the whole time horizon is covered. It consists of the 
following sequence, which is illustrated in figure 

(a) Extend the time horizon by r from to 
Tn+l = Tn + r. 

(b) Extrapolate /opt, old(^) within r, e.g. using 
fnxtr in Mat lab. 
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(c) Adapt the Basis splines to the new time hori¬ 
zon Tn+i and represent /pred,oid(^) in the new 
basis, giving /pred,new(t)- In Matlab one can 
use spap2 to do this. 

(d) Use /pred,new(t) as the initial guess for a 
weighted nonlinear least squares fit over the 
extended time span up to Tn+i- 

(e) Judge the results of the fit based on its re¬ 
duced chi-squared value Xred' below a 

specified bound, continue with an additional 
iteration of steps a)-d), repeating until the full 
region of the data is covered. Otherwise, try 
the following fail-back procedures: 

i. Reduce r, the time by which the time 
horizon is extended, and try again 

ii. Increase the number of Basis splines and 
try again 

iii. Try again using a different starting point. 

If all of those fail, we have to resort to increas¬ 
ing the bound on X^ed* 

3. Post-processing. The following steps are optional 
and were performed manually in cases where we 
desired to improve the fit, or examine its behaviour. 

(a) The optimization over the whole time horizon 
was re-run using different numbers of Basis 
splines for 6{t) and U(t). This served as a 
useful check on the sensitivity of the fit. 

(b) The optimization over the whole time hori¬ 
zon was re-run using a starting point based on 
the previously found optimum plus random¬ 
ized deviations. This examines robustness of 
the final fit. 


E. Error estimation 

To obtain error bars of the time-dependent functions 
we use non-parametric bootstrapping [32]. The process 
is summarized as follow: 



FIG. 7: Parametric bootstrap resampling: Predictions 
for J(t), i(t) and U(t) with error bounds obtained using para¬ 
metric bootstrap resampling, assuming quantum projection 
noise. This can be compared with the error bounds obtained 
from the non-parametric method which are shown in Figure 
l^in the main text. The bounds are tighter for the parametric 
bootstrapping. 


3. Post-process samples 

(a) Form a histogram of the chi-squared values 

^red,r * 

(b) Find and fit a normal-like distribution to the 

histogram with preference to the spread with 
lowest lying x^ed r ^ multi-modal 

distribution. From the fit obtain the mean 
reduced-chi squared value (x?ed r) 

the standard deviation cr^. 

(c) Eliminate the outlier samples by removing all 
dr with Xred r valucs that are 3-5cr^ from the 
mean 

(d) Form a matrix Y where each row vector is 
a sample set of coefficients dr that remained 
after step 3(c). 


1. Estimate initial solution Estimate the time de¬ 
pendent functions from the original data using 
Hamiltonian estimation. 

2. Resampling Create sample solutions for all 
time-dependent functions in the following way: 

(a) Eorm a sample set by randomly picking with 
replacement from the photon count data used 
in qubit detection. 

(b) Re-estimate new time-dependent functions by 
optimizing over the full time span, using the 
solution found in (1) as a starting point. 

(c) Record the reduced chi-squared values Xred r 
for each sample r along with the B-spline 
curve coefficients a^ 


4. Obtain statistics 

(a) Eind the mean B-spline coefficients (d) of 
equation [Ts] by taking the mean over the col¬ 
umn vectors of Y with each element of the 
mean given by {d)i = {ai). 

(b) Eind the covariance matrix X) = coy{Y*^) with 

Y^ij = E [((a^ — ((a^)) {aj — ((aj))] with E the 
expectation operator. The standard devia¬ 
tions of each of the mean coefficients {ai) is 
given by We record these val¬ 

ues in a row vector 

We have also applied parametric bootstrapping in or¬ 
der to obtain the error bounds shown in figure The 
difference to the non-parametric case is that in point (2) 
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a) 


b) 





FIG. 8: Measured data, estimation and residuals: Spin population as a function of detuning and switch-off time of the 
laser beam, for the data sets used to obtain the reconstructed parameters shown in figure]^ a) uses a velocity profile with only 
small variations, b) A second data run in which large variations in the velocity profile were used. Each data point results from 
100 repetitions of the experimental sequence. For the Hamiltonian estimation the data was weighted according to quantum 
projection noise. 


the samples are created using the solutions obtained from 
(1) and adding quantum projection noise. For each sam¬ 
ple the Hamiltonian is estimated. The estimates from 
multiple samples are used to construct error bounds in 
the same manner as for the non-parametric resampling. 
We have found that the error bounds obtained from para¬ 
metric bootstrapping are lower compared to that of the 
non-parametric case as shown in figure We think this 
is due to the latter exploring deviations around a single 
minimum in the optimization landscape, while the case 
resampling arrives at different local minima which are 
spread over a wider region. 


a second pair of data sets in which we take two different 
velocity profiles using the same beam position. This data 
is shown in Figure Also shown are the best-fits ob¬ 
tained from the reconstructed Hamiltonians. The param¬ 
eter variations obtained from the reconstructed Hamilto¬ 
nians for these data sets can be found in the main text in 
Figure]^ The sampling rate of the data in these data sets 
was 2 MHz, resulting in a Nyquist frequency of 1 MHz. 


F. Single beam profile with two different velocity 
profiles. 


As a check that our method is also able to produce con¬ 
sistent results for the Rabi frequency profile, we measured 


















